On the diffusion coefficient of a photon migrating through a turbid medium: a fresh 

look from a broader perspective 
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Does the diffusion coefficient of a photon depend on time t or the probability of absorption fc? To 
find an answer to the question, photon transport in a medium of infinite extent is analyzed using 
the method of moments. It is pointed out that if D is defined so as to make it depend on t or k, it 
will also depend on the experimental conditions; that the parameter k which enters the stationary 
diffusion equation is in general different from that entering the transient version; and that a hitherto 
unused non-Markovian partial differential equation may be used for treating photon transport. 

PACS numbers: 42.25.Dd, 05.60.Cd, 42.30.Wb, 02.70.Uu 



A theoretical basis for quantifying translational diffu- 
sion in terms of a diffusion coefficient (D) was laid in 1905 
by Einstein J3, ; shortly afterwards, he published a sec- 
ond note 0,0], wherein he acknowledged that his formula 
"cannot be applied for any arbitrarily small time" . Ein- 
stein concluded that the diffusion equation (DE), which 
governs the density in coordinate space, as well as the def- 
inition of D that follows from the equation "only holds 
for intervals of time which are large compared with fiB v ; 
the product \xB, called the velocity relaxation time, will 
be denoted here by 1/0 (see below). To go beyond the 
DE it becomes necessary to formulate a so-called trans- 
port equation for describing the density in phase space; 
it is not sufficient, nor indeed advisable (as shown be- 
low), to employ the telegraph equation (TE). The label 
Lorentz-Boltzmann equation (abbreviated as LBE) will 
be used here for a transport equation applied to photons 
and monoenergetic neutrons |4j. Studies based on the 
LBE have given rise to a debate as to whether D de- 
pends on time and/or the probability of absorption by 
traps in the host medium. A recent addition to this ma- 
terial is a series of papers by Alfano, Cai, Lax and Xu 
0- E3- 01 - who concluded that D depends on time but not 
on absorption; they inferred their definition of D with the 
aid of an approximation based on a cumulant expansion 
of the distribution function I(r,c,t). Earlier, Aronson 
and Corngold y| had claimed the converse to be true, 
and cited some experimental evidence in support of their 
claim. My object here is to examine these issues (which 
arise also in other contexts) from a broader perspective. 

It will be convenient to begin by drawing attention to 
some implications of the cumulant approximation which 
appear to have been overlooked by its proponents, for 
whom I will use the alias CLAX (composed of their ini- 
tials). CLAX proposed that the exact particle density 
N(r, t) in a system with plane symmetry (around the z- 
axis) can be replaced by a product of three Gaussians 
A( G )(r, t) = G(x,t)G(y,t)F(z,t). Since 

N (G) 

can be 

recovered by means other than cumulant expansion, I 
will refer to it as the Gaussian approximation (GAP, for 
short). My first task is to derive, using more elementary 
means and allowing for more general initial conditions, 
the expression for F(z, t). 



Let ez and CI be unit vectors along the z-axis and the 
velocity c, respectively, dCl = dfid<p, /i = ez ■ CI, and 
define 



y(z,ii,t) = 2n 



oo poo 



I(r, n,t)dxdy, (1) 



OO «/ — OO 



and 



F(z,t) = 



/OO pOO pi 
/ N(r,t)dxdy= / fi,t)dfi. (2) 

- oo J — OO J — 1 

We will begin by assuming that absorption is absent 
and scattering is isotropic; with a = c/l = f/r, where 
I is the mean free path and r the mean free time, the 
equation of radiative transfer can be written as 

d t + c/xd z ] 9(z, (i, t) = -a [V(z, n, t) - \F{z, t)] . (3) 

Since the speed c is taken to be a constant parameter, 
it will be convenient to call /x the velocity (along the z- 
direction). 

Calculation of moments. We define the expectation 
of a quantity \ as x M0 = /_ x dfif^ dz x^(z, H, t\ Zo, Mo), 

and proceed to calculate z 2 ^° for a photon with an initial 
velocity /xo and initial position zq. If the initial velocity 
distribution differs from S(fi — fj,o), an ensemble average 
(denoted by an overline) can be found by taking a second 
average (over the initial velocities) of the first average: 
X = X Mo ■ This notational device for distinguishing be- 
tween the two averages is borrowed from Uhlenbeck and 
Ornstein (U&O) 9]. For a medium of infinite extent we 
can insist that F(z, t) — > as \z\ — > oo. If we multiply 
Eq. (3) by \ = z rn fi n and integrate over all z and all /z, 
we obtain the equation 



d d f°° f 

dt^° ~ ° = " ar ° + f / dzF ^ 7 1} / ^ X ' 

(4) 

which can be easily solved to get dx^°/dt. For (m,n) = 
(0, f ), (f ,0) and (0,2) one finds 



z»° = z + cfj, (l-e- at 

~ 9^° 2 —at i 1 / i 

u A = fi^e + |(1 - e 



)/«. 

-at\ 



(5) 
(6) 
(7) 
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The last relation is needed when one goes on to find 
and solve, by setting (m,n) = (1, 1), the equation satis- 
fied by z/J Mo . Finally, with (m,n) = (2,0), the required 
equation for z 2 ^° comes out in the form dz 2 ^° jdt = 
2czJI AI0 , which integrates to give 



2t chqZq{1 — e 



-at, , 2 2 r 2 



2 ^ 2 (3 M 2 -1) 



(erf - 1 - 
ate~ at ) . 



(8) 



One can go on and find z nP,a , for n > 2; by availing oneself 
of symbolic computation, one can reduce the tedium and 
the risk of making errors. Once the raw moments are at 
hand, one can find the central moments Z n ^° , where 



Z = z 



(9) 



CLAX have allowed for anisotropic scattering; their ex- 
pressions for the mean and the variance agree with mine 
if one sets g\ = g 2 in the former and \iq = 1 in the latter; 
it will be sufficient to spell out the variance: 



Z**- = l(?T[2t-(l 



Ae- 



rate 



-3e- 2Qt )/a] . 

(10) 

For a Gaussian distribution, all odd central moments 
vanish, and even central moments satisfy the relation 
M 2n = (2n - 1)!![M 2 ]". One can easily verify that 
Z^VOand z4»° ^ 3[^ W ] 2 - 

Introduction of boundaries. For an infinite-medium 
problem, the asymmetry in the distribution can be taken 
into account by calculating the higher moments, and us- 
ing the results to find a more realistic analytical form 
0, However, this strategy cannot be used when 

boundaries are present. Under such circumstances, GAP 
becomes, I suggest, worthy of consideration, because a 
Gaussian of the form 



F{z,t) 



1 



: exp 



{z — m(t)Y 
2a 2 (t) 



satisfies the partial differential equation 

d t F(z, t) = \a{t)d zz F{z, t) - b(t)d z F(z, t) 



(11) 



(12) 



with a = da 2 jdt and b = dm jdt ^J. Of course, there 
still remains the question of inferring the boundary con- 
ditions to be imposed on a solution of Eq. (12); for further 
details, the reader is referred to a recent article [l2| . 

The definition of D. For defining D, it will be helpful 
to enlarge the scope of our discussion and recall the corre- 
sponding expressions for a Brownian particle (B-particle) 
of mass to 





= u exp(-/3t), 




(13) 


z v ° 


= z o + v o (l-e- 0t )//3, 




(14) 




= v 2 e~ 2f3t + {kT/m){l 


-e-^% 


(15) 




= (fcT//)[2t-(3-4e- 


+ e- 2/3t )//3] , 


(16) 


~zJ 


= (kT/f) [2t-2(l-e- 


-*)I0\ . 


(17) 



Here vq = v(0), v is the z-component of the veloc- 
ity (— oo < v < oo), / denotes the friction coeffi- 
cient, [3 = f/m, and the other symbols have their usual 
meanings |9(. U&O obtained these results through the 
Langevin equation, but one can also utilize the Klein- 
Kramers equation fKKE)[l3j and apply the procedure 
outlined above, choosing x = z m v n and replacing inte- 
gration over all fi by integration over all v; whatever the 
route, one finds that the density F(z,t; zq,vo) remains a 
Gaussian at all times. 

The transient terms in Eqs. (5)-(8) and (13)-(17) rep- 
resent, it cannot be overemphasized, the ballistic phase, 
during which the velocity relaxes exponentially to its 
equilibrium value; diffusion, as envisaged by Einstein, 
starts after this relaxation is over. 

CLAX have opted for the following definition: 



dzZ 2 F(z, t; Zq, jtio) = 2Dt. 



(18) 



To be exact, they write 2D zz ct on the right-hand side, 
but I have absorbed the speed c into the definition of D so 
that its dimensions come out to be [L] 2 [T] _1 . I submit 
that their definition is unreasonable because it implies 
that D depends (not only on t but also) on the initial 
conditions. The truth of this assertion, though implicit 
in the definition, will now be made manifest through an 
explicit calculation of the variance for three cases. 

Thus far we have been occupied with particles with the 
same initial velocity (/^o = 1). The CLAX expression for 
D emerges upon dividing the right-hand side of Eq. (10) 
by 2t. Let us now consider two other situations: isotropic 
distribution of initial velocities (7*o = 0, /Xq = 1/3) and 

a 'magic-angle' incidence (/Zq = 1/V3, Mo = 1/3) • in the 
first case, Eq. (6) gives z = zq, so that Eq. (8) implies 



Z 2 



±c 2 

3° 



T[2t-2(l-e- at )/a], (19) 



and in the second 



Z 2 ' 



¥ 2 - 



[2t - (3 



4e" 



-2at 



)/«] 



(20) 



The definition advocated by CLAX would lead to a new 
expression of D in each case. We will see later that the 
standard definition amounts to replacing the left-hand 
side of Eq. (18) by its long-time limit (t — > oo). 

A caution concerning the TE. If one sets \c 2 t = D = 
kT/f and a = /3, Eq. (17) coincides with Eq. (19), and 
with the result obtained by using the TE; this is to be 
expected because the TE can be obtained from the LBE 
or the KKE only if one assumes that the velocity dis- 
tribution corresponds to equilibrium. The shortcomings 
of the TE were well exposed, within the context of the 
KKE, first by Hemmer 0] and later by Wilemski [l5| . 
Those who place their trust in the TE will do well to de- 
vote some attention to these penetrating analyses, which 
are not vitiated by a change in the transport equation. 

Some other special cases. The formal identity of 
Eqs. (16) and (20) is noteworthy, and so is the fact (which 
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can be easily verified) that if one sets (72 = in the CLAX 
result for the variance, it coincides with Eq. (20). Unfor- 
tunately, this coincidence does not extend to the higher 
moments, which means that photons cannot be made to 
mimic Brownian motion by a judicious adjustment of the 
experimental conditions. 

Does D change in an absorbing medium? If the ab- 
sorption probability per unit time is a constant (fc, say), 
the LBE will take the following form (with 7 = a + k): 



d t + cnd z (z, /i, t) = -7#(z, (Jt, t) - haF(z, t). (21) 



When absorption is present, the densities ^(z,/J,,t) and 
F(z, t) cannot be normalized. At first sight, this can be 
easily remedied, since the substitutions 



V(z, fi, t), F{z, t) = e kt F{z, t) (22) 



undo the loss caused by absorption. The claim made 
by CLAX, concerning the insensitivity of D to k is based 
essentially on this transformation. Though ^(z, /U, t) and 
F(z, t) are properly normalized, their introduction is not 
sufficient to clinch the argument. Consider, for example, 
a system where absorption dominates. In this 
particle will be absorbed at its first (or, at most, the 
second or the third) collision, and this is a time domain 
where ballistic behavior dominates the dynamics. I also 
add in passing that the above argument applies equally 
well to the KKE (or any other linear kinetic equation). 

When a B-particle moves through the suspending 
medium, the direction of its velocity does not change 
much after a collision, and many collisions are needed to 
randomize an initially imparted velocity. In photon (or 
neutron) transport, sometimes called inverse Brownian 
motion [Tfil ]. a particle with no or negligible mass col- 
lides with stationary targets, and Jl is randomized after 
each collision (when scattering is isotropic). Nonethe- 
less, if certain restrictions are satisfied, the DE will ap- 
ply to Brownian motion as well as its inverse |l6j . That 
the fundamental solution of the DE is a Gaussian is 
merely an affirmation of the central limit theorem, ac- 
cording to which the probability density function of the 
sum of N independent, identically distributed random 
variables can be approximated, if N is sufficiently large, 
by a Gaussian, regardless of how the individual variables 
are distributed 10]. In other words, one need not worry 
about the nature of the particle and the host medium; 
conversely, the question of what is diffusing in what can- 
not be answered within the confines of the DE. Though 
Einstein was concerned essentially with B-particles, this 
preamble makes it clear that his derivation would remain 
valid, when suitably applied, to any system. 

Let us therefore follow Einstein p|, and assume the 
existence of a time interval, say e, which is very small 
compared to the time interval, say T, over which ob- 
servations are made, but sufficiently long to justify the 
assumption that the displacements suffered by the par- 
ticle in two successive intervals, each of duration e, can 



be viewed as mutually independent. Whence follows the 
equation 

/oo 
F(z- s,t)<p (s,e)ds, (23) 
-00 

in which (f>o(s,e)ds is the probability that a particle will 
go from z — s to z in a time interval e. By virtue of 
the central limit theorem |l0j| , we can ascribe a Gaussian 
form to </>o(s,e) with zero mean and a variance propor- 
tional to e. Next, we expand F(z — s, t) in a Taylor series, 
introduce the notation s k = s k (f>o(s,e)ds, and get 



F{z,t + e) = F(z,t) + ±s 2 d zz F, 



(24) 



which can be immediately transformed to the DE, dtF 
Dd rr F, with 



D 



lim — . 

e^o 2e 



(25) 



The limit e — > is purely formal, since we have already 
agreed that £>t (or as ^> 1). In fact, it proves advanta- 
geous to replace the limit e — > with e — > 00, for D then 
turns out to be the integral (t = to t = 00) of the ve- 
locity autocorrelation function one gets D = c 2 r/3 
(photons) and D = kT/f (B-particles). 

Let us now apply the same reasoning to an absorbing 
medium (or to a non-absorbing medium but with a diffus- 
ing particle that is subject to first-order decay), and see 
if we can infer a modified DE by replacing (f>o(s, e) with 
cj>k{s,e), the new transition probability. To make any 
progress, it seems necessary to assume that 4>k = <fioe~ 
and then we obtain 



-fee 



F(z,t + e)-F(z,t) = -{\-e- ke ) + \s 2 e- ke d zz F(z,t). 

On dividing by e and letting e — > 0, we formally get 

d t F(z, t) = -kF(z, t) + DO zz F(z, t), (27) 

but the limit is to be interpreted with some sensitivity, 
since we cannot transgress the condition > 1. The 
difficulty can be avoided — and, as we shall see, meaning- 
ful experiments made possible — if we impose the demand 
before taking the limit. We conclude, then, that 
Eq. (27), with the customary definition of D |17|. is valid 
only if k <C a (i.e., negligible absorption). 

I elaborate on the above argument by using an analogy. 
Time-resolved photo-induced anisotropy measurements, 
which monitor the rotational diffusion of an electronically 
excited probe, are interpreted through the rotational DE. 
Though the definition of anisotropy incorporates a nor- 
malization akin to that in Eq. (22), the experiment be- 
comes inconclusive, owing to a precipitous decline in the 
signal-to-noise ratio, if fc 3> T> r , where the rotational dif- 
fusion coefficient T> r is the analog of a. The lifetime of 
the probe 1/k plays the role of T, effectively setting a 
scale for the observation time; the only remedy, when 
diffusion is slow, is to use a long-lived probe [lSl Il9l|. 
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Time- dependent and stationary versions of the DE. 
Aronson and Corngold Q argue that D depends on k, 
and base their claim on analyses of stationary versions of 
the LBE. However, these analyses lead not to D as such, 
but to a quantity, called the diffusion length L, which 
is a measure of the average distance a particle would 
travel, when released in an infinite medium, before it is 
captured. Evidently, L would decrease with increasing 
absorption, and if one chooses to define D through the 
relation D = kL 2 , D too will depend on k. But a D so 
inferred would pertain, if absorption is high, to a parti- 
cle that is captured before it has taken a large number 
of diffusive steps; in the extreme case k ^> a the par- 
ticle would be absorbed during its ballistic phase, and 
such a definition would be liable to an objection similar 
to that leveled against a time-dependent D, since the fi- 
nal form of the Poisson type equation depends also on 
the distribution of the sources [2£| • "It is also necessary 
to remark," Davison [2fJ informed his readers and I go 
along with him, "that, though the value of c [our a/7] 
in the system does not enter directly the criteria of ap- 
plicability of the diffusion approximation, in practice the 
dimensions of the system will usually be of the order of 
L . . . , and this effectively limits the application of the 
theory to systems where |1 — c| is small". 

The time-independent DE corresponds to a stationary 
(or steady) state, a situation when a constant rate of pro- 
duction of the particles within a region equals their rate 
of absorption. When this condition is fulfilled, the time- 
dependent DE, d t F(z,t) = Dd zz F(z,t)-kF(z,t) + S(z), 
where S(z) denotes the source term, goes into the sta- 
tionary version Dd zz F s (z) — kF s (z) = —S(z), which 
leads, upon division throughout by D, to a Poisson type 
equation, formally identical with that deducible from the 



stationary LBE: d zz F s (z) — F s (z)/Lq = —s(z), where 
Lq = D/k and s(z) = S(z)/D. Smoluchowski was the 
first to realize that if a particle is produced at a random 
point in a medium containing traps, the trapping _proba- 
bility per unit time is, in general, a function of t |l3tl21| . 
In other words, absorption or trapping differs from first- 
order decay (such as electronic de-excitation or nuclear 
disintegration). Aronson and Corngold [8| did not take 
this complication into account when they wrote: "It is 
obvious (and is also well known) that the solution of 
the time-independent diffusion equation is obtained from 
that of the time-dependent equation by integrating over 
all time." It is far from straightforward — unless absorp- 
tion is small — to relate the parameter k appearing in the 
stationary DE to the corresponding time-dependent pa- 
rameter which should appear in the time-dependent DE. 
The reader is referred to two reviews 0,H2 and a recent 
article |23| . 

Conclusions. The principal conclusions of this study 
will now be stated:- (a) Since the DE implies and is im- 
plied by the central limit theorem, an unequivocal defi- 
nition of D is possible only if the diffusing particle can 
execute a large number of diffusive steps; when this con- 
dition is satisfied, D becomes independent of time t and 
absorption probability k; if this restriction is ignored, D 
will depend not only on t and k, but also on the ex- 
perimental conditions, (b) The Gaussian approximation 
implies that photon transport can be described by a non- 
Markovian partial differential equation, (c) Inferences 
drawn from the TE are suspect, (d) Smoluchowski's ob- 
servation that the absorption probability per unit time 
depends on t should be taken into account when the DE 
is compared with its stationary version. 
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